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ABSTRACT 

We introduce a new one-dimensional stellar evolution code, based on the existing Dartmouth code, that 
self-consistently accounts for the presence of a globally pervasive magnetic field. The methods involved in 
perturbing the equations of stellar structure, the equation of state, and the mixing-length theory of convection 
are presented and discussed. As a first test of the code's viability, stellar evolution models are computed 
for the components of a solar-type, detached eclipsing binary (DEB) system, EF Aquarii, shown to exhibit 
large disagreements with stellar models. The addition of the magnetic perturbation corrects the radius and 
effective temperature discrepancies observed in EF Aquarii. Furthermore, the required magnetic field strength 
at the model photosphere is within a factor of two of the magnetic field strengths estimated from the stellar 
X-ray luminosities measured by ROSAT and those predicted from Ca II K line core emission. These models 
provide firm evidence that the suppression of thermal convection arising from the presence of a magnetic field 
is sufficient to significantly alter the structure of solar-type stars, producing noticeably inflated radii and cooler 
effective temperatures. The inclusion of magnetic effects within a stellar evolution model has a wide range of 
applications, from DEBs and exoplanet host stars to the donor stars of cataclysmic variables. 
Subject headings: binaries: eclipsing — stars: evolution — stars: individual (EF Aquarii) — stars: interiors — 
stars: low-mass — stars: magnetic field 



1. INTRODUCTION 

The vast array of physics incorporated in standard low- 
mass stellar evolution models (see e.g., Chabrier & Baraffe 
1997; Baraffe et al. 1998; Dotteretal. 2007, 2008, and ref- 
erences therein) appears to be insufficient for predicting the 
properties of low-mass stars. Studies of detached eclipsing 
binary (DEB) systems allow for a very precise determination 
of the mass and radius of the individual stellar components 
with uncertainties commonly below 2%. Over the past two 
decades, these precision studies have accumulated strong ev- 
idence that the radii predicted by low-mass stellar evolution 
models are deflated compared to the observations, at a given 
mass (e.g., Popper 1997; Torres & Ribas 2002; Ribas 2006; 
Morales et al. 2009; Torres et al. 2010; Kraus et al. 2011). 
Typically quoted is that DEB radii appear to be systemati- 
cally 10% larger than model predictions and DEB effective 
temperatures are 5% cooler than the models imply. Further 
evidence has been garnered by studies of single low-mass 
stars (Berger et al. 2006; Morales et al. 2008), which confirm 
the aforementioned trends. Although, Boyajian et al. (2012) 
present results that may be interpreted as counter to these 
claims. 

Recent work has shown that the disagreements may not be 
as severe as previously believed, when the age and metallic- 
ity of the DEBs are considered (Feiden & Chaboyer 2012). 
Still, many DEB systems display moderately inflated radii 
(less than 5%) with a small subset displaying radically in- 
flated radii (upward of 10%) compared to stellar models 
(Feiden & Chaboyer 2012; Terrien et al. 2012). This modest 
radius offset between observations and models and the pres- 
ence of strongly inflated stars, suggests that stellar evolution 
models must invoke new physics to account for the appear- 
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ance of inflated radii. 

Implicated as the culprit for the observed inflated radii 
and cooler effective temperatures is magnetic activity (Ribas 
2006; Lopez-Morales 2007; Morales et al. 2008). The sys- 
tems at the heart of the problem are often close binaries with 
short orbital periods. Tidal synchronization of the compo- 
nents acts to spin-up the rotation rate of each star, enhanc- 
ing the dynamo mechanism and thus supporting a stronger 
magnetic field within each. While there are a number of di- 
rect magnetic field measurements for single stars, there are 
few among fast rotating binaries. Instead, the hypothesis is 
reinforced by the presence of strong chromospheric Ha emis- 
sion (Morales et al. 2008; Stassun et al. 2012), chromospheric 
Ca II H and K emission (Skumanichet al. 1975), as well as 
coronal X-ray emission (Lopez-Morales 2007) among inflated 
stars. The presence of such emission features is considered to 
be the result of the dissipation of magnetic energy in the stel- 
lar atmosphere. 

There are, however, stars from DEBs that display inflated 
radii, despite existing in long-period systems. Both LSPM 
J1112+7626 (Irwin etal. 2011) and Kepler-16 (Doyle et al. 
2011; Winn et al. 2011; Bender et al. 2012) have orbital pe- 
riods of about 41 days, suggesting that the components are 
tidally unaffected by the presence of their companion. It is 
possible that these stars have not had sufficient time to shed 
angular momentum (Skumanich 1972), preserving strong 
magnetic fields that they likely possessed near the zero-age 
main sequence. Whether the low-mass stars in these systems 
are still magnetically active enough to affect the stellar radius 
remains unclear. The work by Winn et al. (2011) appears to 
suggest that Kepler-16 is still relatively young (2-4 Gyr) 
as the primary is more active than the Sun, based on Ca II 
line emission. This supports the notion that magnetic fields 
may be influencing the structure of each component. On the 
other hand, the age of LSPM Jl 112+7626, estimated from 
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gyrochronology, reveals that it is approximately 9 Gyr old 
(Barnes 2010), supporting the notion that magnetic activity is 
likely not playing a prominent role in either star's evolution. 

Further complicating our picture of low-mass DEBs, KOI- 
126 (Carter et al. 2011) curiously matches standard stel- 
lar evolution models (Feiden et al. 201 1; Spada & Demarque 
2012). The mass of KOI-126 B and C and their orbital period 
are very similar to CM Draconis, one of the quintessential 
low-mass DEB systems displaying inflated radii (Lacy 1977; 
Morales et al. 2009; Terrien et al. 2012). It is interesting, 
then, that stellar evolution models would even come close to 
accurately predicting the observed stellar properties of KOI- 
126. Given what is known about CM Draconis, one would 
expect KOI-126 to display inflated radii as a consequence of 
moderate magnetic activity. Although, MacDonald & Mullan 
(2012) speculate that KOI-126 should, in fact, not be terribly 
active and find it unsurprising that standard models should fit 
the system so well. 

Despite the identification of a potential culprit responsible 
for inflating the radii of DEB components, only ad hoc proce- 
dures for treating the effects of magnetic fields have been in- 
troduced (Mullan & MacDonald 2001; Chabrier et al. 2007). 
The method examined by Chabrier et al. (2007) included ar- 
tificially decreasing the convective mixing-length parameter, 
so as to mimic the effect of a global magnetic field within 
the star, as well as artificially reducing the star's bolomet- 
ric flux in an effort to reproduce the effects of photospheric 
spots. A second method, proposed by Mullan & MacDonald 
(2001), altered the Schwarzschild criterion by perturbing the 
adiabatic gradient in a manner consistent with the work of 
Gough & Tayler (1966). 

Investigations by both groups appear to be at odds 
with one another. Chabrier et al. (2007) and Morales et al. 
(2010) claim that starspots appear to be the dominate 
mechanism inflating stellar radii, and that modifications 
to convection require unrealistic magnetic field strengths 
(i.e., reductions in the mixing length in their formula- 
tion). On the other hand, Mullan & MacDonald (2001) and 
MacDonald & Mullan (2012) conclude the opposite that re- 
duction in convective efficiency is ultimately the dominant 
mechanism. Regardless of which is really the dominant mech- 
anism, both approaches are inherently ad hoc, yet both are 
capable of reproducing the observed inflated stellar radii. 

In this paper, we introduce a self-consistent treatment of 
a globally pervasive magnetic field embedded in the frame- 
work of the Dartmouth stellar evolution code (Dotter et al. 
2007, 2008). Our approach follows the outline provided by 
Lydon & Sofia (1995), though we deviate from their method 
in a number of ways that are described below. All of the 
stellar structure equations, including those in the equation of 
state, are self-consistently modified, as opposed to arbitrarily 
altering a single quantity. In this way, modifications to the 
efficiency of thermal convection are accounted for in a more 
complete fashion, owing to the full thermodynamic treatment 
of the magnetic field. Overall, the approach used to model 
magnetic effects can be considered analogous to the parame- 
terized mixing-length treatment of convection. 2 The viability 
of the models is tested against results from a recent study that 
characterized the DEB EF Aquarii (Vos et al. 2012). 

EF Aquarii (HD 217512; henceforth EF Aqr) is a solar- type 
DEB found to contain two components displaying drastically 

2 In so far as reducing an inherently nonlinear, three-dimensional process 
into terms suitable for a one-dimensional model. 



Table 1 

Fundamental Stellar Parameters for EF Aqr 



Parameter 


EF Aqr A 


EF Aqr B 


M(M ) 


1.244±0.008 


0.946±0.006 


R(R @ ) 


1.338±0.012 


0.956±0.012 


Tiff (K) 


6150±65 


5185±110 


[Fe/H] 


0.00±0.10 



inflated radii (Vos et al. 2012). Fundamental parameters of 
the system are quoted in Table 1. What is most striking, is the 
similarity of the secondary to the Sun and the entire system to 
a Centauri A and B, in terms of the stellar masses and compo- 
sition. Although the secondary appears similar to both a Cen 
B and the Sun, its radius appears to be about 10% larger than 
one would expect based on stellar evolutionary calculations. 
The effective temperature of the primary further reveals that 
both components suffer from substantial radius inflation. 

In an effort to reconcile the observations with predictions 
from theoretical models, Vos et al. (2012) reduced the value 
of the convective mixing-length. They found OCmlt of 1.30 
and 1 .05 were required for the primary and secondary, respec- 
tively, compared to their solar-calibrated value of 1.68. They 
concluded that fine-tuning the models allows for an accurate 
description of the observed properties. 

Reduction of the required convective mixing length may 
be physically motivated in two ways: (1) naturally ineffi- 
cient convection and (2) magnetically suppressed convection. 
While we must be careful to not read too much into the re- 
ality of mixing-length theory, in stellar evolution models the 
mixing length is an intrinsic "property" of convection. Thus, 
reducing the mixing length is akin to saying convection is not 
very efficient at transporting excess energy. 

Bonacaetal. (2012) calibrated the convective mixing 
length for solar-like stars using asteroseismic results provided 
by the Kepler Space Telescope. They found that the value of 
OCmlt in stellar models is tied to stellar properties (i.e., \ogg, 
logr e ff, and [M/H]). Applying the Bonacaetal. (2012) em- 
pirical calibration to the stars in EF Aqr, we find that the 
primary and secondary component require OCmlt = 1 -68 and 
1.44, respectively. Again, compared with their solar cali- 
brated value of OCmlt = 1-68. The asteroseismically adjusted 
mixing-lengths are significantly larger than the fine-tuned val- 
ues determined by Vos et al. (2012). Therefore, it appears that 
naturally inefficient convection is insufficient to explain the 
inflated radii of EF Aqr. We are left with the option that mag- 
netic fields may possibly be to blame. 

In what is to follow, we describe a self-consistent approach 
to modeling the effects of a globally pervasive magnetic field 
with application to the EF Aqr system. Details of the stan- 
dard Dartmouth models are presented in Section 2 followed 
by a description of the magnetic perturbations introduced to 
the code in Sections 3 and 4. Section 5 demonstrates the abil- 
ity of the invoked perturbations to reconcile the models with 
the observations. We conclude with a further discussion of the 
results and their implications in Section 6. 

2. MODELS 

The framework throughout which magnetic effects are 
invoked is provided by the Dartmouth Stellar Evolution 
Program (DSEP; Dotter et al. 2008), 3 a descendant of 

3 Available at: http://stellar.dartmouth.edu/models/ 
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the Yale Rotating Evolution Code (Guenther et al. 1992). 
Standard stellar evolution henceforth refers to the ba- 
sic physics without any magnetic perturbation. The 
standard physics incorporated in DSEP has been de- 
scribed previously (Chaboyer & Kim 1995; Chaboyer et al. 
2001; Bjork & Chaboyer 2006; Dotteretal. 2007, 2008; 
Feiden et al. 2011), although we will briefly review the ele- 
ments that are pertinent for the current study. 

Above 0.80M G , DSEP invokes an ideal gas equation of 
state (EOS) supplemented by a Debye-Hiickel correction in 
order to account for ion-charge shielding (Chaboyer & Kim 
1995). This EOS is computed analytically, in a self-consistent 
manner, within the code and does not rely on the interpo- 
lation within EOS tables. Opacities are drawn from two 
sources, the OPAL opacities for the high temperature regime 
(Iglesias & Rogers 1996) complimented by the Ferguson low- 
temperature opacities (Ferguson et al. 2005). Surface bound- 
ary conditions are defined using the PHOENIX AMES-COND 
model atmospheres (Hauschildt et al. 1999a,b), attached to 
the model interior where T = T e g. 

Atomic diffusion and the gravitational settling of helium 
and heavy elements are implemented using the prescription 
of Thoul et al. (1994). Additional diffusion effects associated 
with turbulent mixing (Richard et al. 2005) are also included. 
Details of the latter are presented in Feiden et al. (201 1). Fi- 
nally, convective core overshoot is treated following the meth- 
ods outlined by Demarque et al. (2004). 

Required before any analysis pertaining to stellar evolution 
models is calibrating the model properties to the Sun. Deter- 
mination of the initial solar helium and heavy element mass 
fractions along with a compatible mixing-length parameter 
(Finit, Zi n it, and ocmlt, respectively) was performed by cali- 
brating a 1 Mr., model to the Sun. Our solar model was re- 
quired to reproduce the solar radius, solar luminosity, radius 
at the base of the convection zone, and the solar photospheric 
(Z/X) at the solar age (4.57 Gyr; Bahcall et al. 2005). The 
final set of parameters necessary to satisfy the above crite- 
ria for the Grevesse & Sauval (1998) solar composition was 
Finn = 0.27491, Zi n i t = 0.01884, and oc M lt = 1-938. 

3. MAGNETIC PERTURBATION 

3.1. Magnetic Field Characterization 

Investigating the effects of a global magnetic field on the 
interior structure of a star over long time baselines, requires 
formulating a purely three-dimensional (3D) phenomenon in 
terms suitable for a one-dimensional (ID) numerical model. 
Unfortunately, full 3D magnetohydrodynamic (MHD) models 
are not yet capable of modeling stellar magnetic fields over 
the long time baselines required for stellar evolutionary cal- 
culations. This is in part due to the rapid 4 diffusion of the 
magnetic field and the immense computational time required. 
Therefore, in order to probe the effects of a magnetic field, we 
seek to avoid directly solving the induction equation 



^ = Vx(uxB)- 
dt 



r|V 2 B. 



(1) 



While not actively seeking a solution to the full suite of 3D 
MHD equations, it is possible to use the theoretical frame- 
work of MHD to provide a reasonably accurate ID descrip- 
tion of a magnetic field and its associated properties. Ulti- 
mately, we are able to describe a magnetic field in terms of 

4 Relative to a typical stellar lifetime. 



the MHD equations and then project out the radial compo- 
nent, the component necessary for stellar evolutionary model 
computations. 

The spatial and temporal evolution of a given magnetic field 
are governed, quite naturally, by Maxwell's equations, 



V-E = 

1 8B 

VxE = — ^~ 
c dt 

V-B = 
„ 471 

V x B = — J 

c 



(2) 
(3) 
(4) 
(5) 



where within the stellar plasma, we assume any regions of 
excess charge inducing an electric potential will rapidly neu- 
tralize owing to the mobility of other charges (Debye shield- 
ing). Thus, we can safely assume that the plasma is electri- 
cally neutral, p e = 0. For simplicity, we here made another 
assumption, that temporal variations of the large-scale field 
are small, suggesting that the conduction current dominates 
the displacement current. 

Now, let us consider the interactions between the elec- 
tric and magnetic fields within a dense, ionized fluid mov- 
ing with arbitrary velocity, u. For slow temporal evolution, 
non-relativistic dynamics may be described by a single con- 
ducting fluid that obeys the classical equations of hydrody- 
namics coupled with the equations of electromagnetism; the 
MHD equations (Jackson 1998). Considering a perfectly con- 
ducting, non-viscous, non-rotating, compressible fluid in the 
presence of a gravitational field, the MHD equations govern- 
ing the system are Ohm's law for a moving fluid, 



J = o E 



uxB 



the equation of mass continuity. 

dp 



d, 



V-(p m u)=0, 



(6) 



(7) 



where p,„ is the mass density, and the fluid equation of motion, 



du JxB 



Pm dt 



-V.V- 



Pmg 



(8) 



with g being the gravitational field vector and P representing 
the gas pressure tensor. The electromagnetic term in the fluid 
equation of motion is associated with the assumption that a 
magnetic field permeates the plasma. However, we have ne- 
glected forces associated with any electric fields, for reasons 
detailed above. 

Since, a priori, we have no knowledge of the current den- 
sity within a given fluid, we replace the current density within 
Equation (8) using Equation (5). The equation of motion may 
now be written as 



Pm^r = j-(VxB)xB-V.V + p m g. 
dt 4tc 



(9) 



With the aid of a vector operation identity and knowing that 
the magnetic field is divergenceless, this may again be rewrit- 
ten as 



p m ^ = ±(B.V)B--VB 2 -V.V- 
Vm dt 4tc v ; 8tc 



Pmg- 



(10) 



Immediately, we recognize that the electromagnetic contribu- 
tions on the right-hand side are the familiar magnetic tension 



4 



Feiden & Chaboyer 



and pressure terms, respectively. However, the final form of 
the equation of motion requires one further step. The first 
two terms on the right-hand side may be expressed as the di- 
vergence of a magnetic stress tensor (Gurnett & Bhattacharjee 

2005), let us call it V, such that 



fc4 BB <— >fi 2 

t =-~r+ i 



(11) 



with I representing the identity tensor. This definition then 
implies, 



Ping- 



(12) 



The magnetic stress tensor introduced above can be thought 
of an anisotropic pressure tensor, where the pressures it de- 
scribes are intrinsic properties of the magnetic field. 

Stars, however, are considered to be in hydrostatic equilib- 
rium. This implies the absence of bulk fluid motion, forcing 
the left-hand side of the equation of motion to vanish. There- 
fore, 



= Ping-, 



(13) 



which is a statement of magnetohydrostatic equilibrium. 
3.2. Stellar Structure Perturbations 

At the most fundamental level, one-dimensional stellar evo- 
lution codes simultaneously solve a set of four coupled, first- 
order differential equations. 5 They are the equation of mass 
conservation, hydrostatic equilibrium, energy transport, and 
energy conservation. Qualitatively, we can easily predict how 
these equations will be altered by the presence of a magnetic 
field which may then be translated into a quantitative descrip- 
tion. 

The equation of mass conservation should be unaltered by 
any magnetic perturbation. Of course, this is assuming that 
mass removed by stellar winds is negligible and that transient 
events that may remove mass (i.e., flares, coronal mass ejec- 
tions) are neglected. The stated conditions hold for our ap- 
proach. Thus, 

dr 1 

— = (14) 

dm 4%r 2 p 

where we have dropped the subscript m on the density and 
assume all references to density are specifically to the mass 
density, unless otherwise noted. 

Hydrostatic equilibrium, as we saw earlier in Equation (13), 
is modified through the inclusion of the magnetic pressure and 
tension. Projecting out the radial component of the magnetic 
pressures, we are able to adapt the three-dimensional concept 
for one-dimensional models. Therefore, we have 



dP 
dm 



1 



Gm 



4nr 4nr z p 



(B • V)B 
4k 



(15) 



The precise handling of the vector magnetic field within the 
code will be discussed later. 

The final form of the energy transport equation is the same 
as if there were no perturbation. Namely, 



dT 
dm 



= -V 



temp 



dp 

dm 



(16) 



5 There are additional equations often included to account for atomic dif- 
fusion. While included in DSEP, we do not seek perturbations to these equa- 
tions at the present time. See Mathis & Zahn (2005) for a rigorous treatment 
of mixing associated with magnetic fields. 



where V temp is the local temperature gradient. Magnetic 
perturbations to the stellar structure equations will self- 
consistently alter the temperature gradient through various 
thermodynamic considerations. Of greatest importance will 
be the affects on the treatment of convection. The full treat- 
ment will be discussed in the next subsection. 

Finally, there are changes to the parameters present in the 
canonical equation of energy conservation in stellar evolution. 
Modifications to these parameters arise from the treatment of 
the specific thermodynamic equations (discussed in the next 
subsection) and additional terms that are electromagnetic in 
origin. The final form of the energy conservation equation is 



dL 

dm 



dU 
~dt 



P_dp t Q, 
p 2 dt 



ohm 



fpoynt 



(17) 



Aside from the first three standard terms on the right-hand 
side, there are two additional electromagnetic terms. First, 
there is a Poynting flux associated with the field, 



oynt 



— ExB, 

471 



(18) 



although as discussed above, we assume the Zs-field is zero 
everywhere. Next, energy is also associated with the Ohmic 
dissipation of electric currents brought about by the resistive 
nature of the plasma. 



ohm 1 



I 2 R, 



(19) 



where / is the electric current and R is the resistance of the 
medium. Here, electrical currents are converted to heat that 
then is transmitted to the surrounding plasma. Since we have 
assumed an infinitely conducting plasma, this energy term 
goes immediately to zero. 

3.3. Thermodynamic Considerations 

The effects of a global magnetic field are introduced into 
the thermodynamic framework supplied by DSEP follow- 
ing the approach outlined by Lydon & Sofia (1995, hereafter 
LS95). Providing a detailed, step-by-step guide of the mag- 
netic perturbation to the various thermodynamic quantities 
would prove tedious. Therefore, we refer the reader to LS95 
for a full derivation of each equation presented below. Our 
aim in this subsection is to adequately summarize the perti- 
nent aspects of LS95 and highlight where we diverge from 
their original approach. 

At the core of the LS95 method is the specification of a new 
thermodynamic state variable, %, such that 



(20) 



The state variable % is the magnetic energy per unit mass and 
B{r) is the magnetic field strength at radius r. Unlike LS95, 
our definition of % depends on the radial distribution of the 
magnetic field strength and also on the density of the stel- 
lar plasma. Originally, LS95 favored a mass-depth-dependent 
function, %(M r ). However, we moved away from this pre- 
scription when we realized several thermodynamic derivatives 
became divergent. Once the magnetic field strength is speci- 
fied throughout the star, it is straightforward to calculate % at 
each point within the model. 

The energy associated with the magnetic field arises due to 
forces exerted by the magnetic field on the plasma. These 
forces are represented by the anisotropic pressure tensor 
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present in Equation (12). As a first approximation, we con- 
vert the pressure tensor to a scalar pressure by taking the trace 
of the pressure tensor to yield the mean magnetic pressure, 



1 



BB <->B 2 



(21) 



Since we are not solving the full set of MHD equations, we 
look, instead, to set approximate upper and lower limits on the 
scalar pressure. Assuming a Cartesian coordinate system, if 
we imagine the magnetic field is parallel to the z-axis, or for 
a star, the rotational axis, then we may expand the pressure 
tensor to read 



V 



B 2 /8n 
B 2 /Sn 







-B 2 /4% + B 2 /8n 



(22) 



Note that there is an isotropic magnetic pressure associated 
with each diagonal element along with the additional mag- 
netic tension term in the final element. Since tension is di- 
rected along the field line, the tension exists in the z-direction 
only, in this instance. Taking the trace, we find 



(-Pmag) 



3 \d,n) 



1 



:XP- 



(23) 



The above equation is satisfied for a magnetic field where 
a strong tension component is present. However, if we as- 
sume that there is no tension at all, then, following the same 
procedure, 

^a g ) = 3( W j=Xp. (24) 

We can now limit the strength of the scalar magnetic pressure 
within the one-dimensional framework. Specifically, 



1 

3XP < (J a mag)<ZP- 



(25) 



Defining a "geometry parameter," akin to LS95, allows us to 
emulate the effects of having a strongly curved field or a field 
with no curvature, and varying degrees between the two ex- 
tremes. This geometry parameter is defined such that we re- 
cover the average magnetic pressure for each case above, 



where 



_ ( 2 tension-free 
^ 1 4/3 maximum tension 



(26) 
(27) 



In both cases, the appropriate expression for the scalar mag- 
netic pressure is returned. With the magnetic energy density 
and pressure formulated as scalars, we have successfully con- 
verted the inherently three-dimensional magnetic field into a 
one-dimensional magnetic perturbation. In the process, we 
have also reproduced the scalar parameters originally pre- 
sented by LS95. 

3.3.1. Equation of State 

The derivations that follow hereafter in Sections 3.3.1 and 
3.3.2 are provided as a review of the LS95 method to enable 
transparency and enhance the clarity of discussions concern- 
ing the application of our models. Original, complete deriva- 
tions are to be found in LS95. We do, however, deviate from 



their paper in Equation (65), where it is stated explicitly be- 
low. 

We have just seen that magnetic fields exert forces on the 
plasma and, thus, carry an associated pressure, tension, and 
energy. The introduction of these terms into the equations of 
stellar structure then necessitates the inclusion of these pa- 
rameters in the EOS of the system. Again, we will mention 
only the most important modifications, deferring to LS95 for 
a rigorous treatment. Beginning with the first law of thermo- 
dynamics, 

dQ = TdS = dU + PdV (28) 

we recognize that each term contains, now, both the standard 
gas and radiation terms as well as a new magnetic contribu- 
tion, 

dQ = T(dS + dSy) = {dU + dU x ) + P dV. (29) 

In the above equation, the magnetic perturbation rightly does 
not contribute any work. However, in order to write the equa- 
tion as a function of the total pressure, we may subtract off 
the magnetic contribution, 

TdS = dU + PdV -(y-l)^dV (30) 

where we take the volume to be the specific volume, V = p _1 . 
Hereafter, it is also assumed that any unsubscripted quantity 
refers to the total quantity while gas and radiation are lumped 
under the subscript (zero) convention and magnetic vari- 
ables carry a subscript %. 

Equation (30) is the new "non-standard" first law of ther- 
modynamics and suggests that 



and 



Vorp=f(P,T, X ) 



UorS = f(p,T, x). 



Following the derivation by LS95, explicitly writing out the 
other state variables illustrates the effects of adding the mag- 
netic perturbation. Ignoring constants for clarity and ease, 



p = P r + ir 4 +xp(y-i) 

p=[P-T 4 /3]/[T + (y-l) 



U- 



3 T 4 
2 p 



which are all subject to the EOS, 

dp dP <,dT d% 

— = a v— . 

p P T x 

The coefficients in the EOS above are defined as follows, 

{d\nPj Tx 

§ = /ainp\ 

\d\nTj Rx 

V = -(^ 

\d\nxJ PT 



(31) 
(32) 

(33) 

(34) 

(35) 
(36) 
(37) 



and will be referred to, as such, throughout the rest of the 
paper. Note, that a carries no subscript and should not be 
confused with the convective mixing-length parameter, OCmlt- 
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An immediate consequence of altering the thermodynamic 
variables is the effect on the specific heats, 



c P - 



c v = 



dQ 
dT 

dQ 
dT 



P,X 



dU 
dT 

dU 
dT 



v,x v"- 1 / v,% 
which are related to one another via the relation, 



Cp-Cy = 



PS 2 
pTa' 



(39) 



(40) 



The difference in specific heats is written just as it would 
be without any magnetic perturbation, however, each term is 
self-consistently modified by the presence of a magnetic per- 
turbation. 

The change in heat as a result of the magnetic perturbation 
follows and is found to be 



dQ = c P dT - -dP 



/P8v 



1 Ux, 



(41) 



although, the addition of the purely magnetic term, d%, as a 
result of the perturbation should not be included. The reason- 
ing for this is simple. If we assume that magnetic phenomena 
are generated through the dynamo action, then we inference 
that rotational energy is the source for the energy converted 
to the pure magnetic term. Since our models do not account 
for rotation, we discard the final term in the above equation. 
Thus, the change in heat to be considered in the stellar lumi- 
nosity equation is 



8 P8v 

dQ = c P dT - -dP H d%. 

P P«X 



(42) 



LS95 were quick to point out, that at the instant of any mag- 
netic perturbation, the change in heat due to magnetic effects 
should be exactly zero. 

Finally, from the observed change in heat comes the defini- 
tion of the adiabatic gradient. Adiabaticity requires constant 
entropy, and therefore no heat exchange, meaning 



P8 



'ad 



/dlnT\ 
" \ dlnPj Sx ~ pTcp 



(43) 



as it is in the non-magnetic case. Again, only in appearance; 
the actual variables are altered by the introduction of a mag- 
netic perturbation. 

3.3.2. Mixing-length Theory 

Convection is determined to occur in regions where a given 
fluid parcel is unstable to a small displacement in the radial 
direction. The primary method of determining convective sta- 
bility is to analyze the density of a generalized fluid parcel. 
Parcels that are less dense than their surroundings will travel 
radially outward until they reach a height within the star at 
which the surrounding fluid has the same density as the par- 
cel. 

Upon reaching this point, the fluid parcel is assumed to fully 
mix with its surroundings becoming indistinguishable from 
the rest of the fluid. Conversely, if a fluid element is more 
dense than surrounding fluid, it will sink down to a greater 
depth in the star, following the same trend as a rising convec- 
tive element. In either case, gravity is the restoring force. One 



assumption is that the fluid parcel is considered to always be 
in pressure equilibrium with its surroundings. 

The distance over which a fluid parcel travels before mix- 
ing is the well-known "mixing-length." Mixing-length theory 
(MLT) has been well established as a local means of prescrib- 
ing convection for a one-dimensional stellar evolution code 
(Vitense 1953; Bohm-Vitense 1958). At locations where vari- 
ous differences in prescriptions of MLT occur, we will specify 
our assumptions. 

Stability of a fluid parcel is determined by comparing the 
density of the element to that of the surroundings 



Dp = 



Ar 



(44) 



where e and s denote quantities of the fluid element under 
consideration and the surroundings, respectively. To ensure 
stability, we require Dp > 0, meaning, the element is stable to 
small radial displacements, 



dp 
dr 



>0. 



(45) 



As the fluid parcel is displaced radially, LS95 reminds us that 
we must consider how % of the parcel reacts. If the initial % of 
the parcel does not change as the element is displaced, 



( d\ni\ 



V dr 



= 0. 



Conversely, if % of the element is always equal to that of the 
surrounding material, there must be a flux of % as the element 
is displaced, 

{ d\n%\ = ( d\n% \ 
~L \ dr : 



V dr 



It is therefore advantageous to relate the spatial gradient of 
magnetic energy density of the parcel to that of the surround- 
ings by introducing a free parameter, /, such that 



d\n%\ 
dr 



(46) 



where / varies between and 1. Later, we will attempt to 
eliminate this free variable and set it to a physically realistic 
value. 

Expanding the density differentials introduced in Equa- 
tion (45) and multiplying through by the pressure scale 
height 6 casts the stability criterion according to known dif- 
ferentials, 



8V,-(l-/)vV z -8V tem p>0. 
where we made use of a series gradient definitions, 



^temp — 



/d\nT\ 
fdlnT\ 

\d^p) e 

_fd\n X \ 
\d\nP J s 



(d\n X \ f dr \ 
V dr JAd\nPj; 



(47) 



(48) 



(49) 



(50) 



Here, we reveal that we are basing our MLT formulation on the pres- 
sure scale height, Hp = —(dr/dlnP), as opposed to the temperature scale 
height since we assume pressure equilibrium between the fluid parcel and its 
surroundings. 
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Here, we have defined the temperature gradient of the sur- 
rounding plasma, temperature gradient fluid element in ques- 
tion, and the magnetic energy density gradient, respectively. 

The magnetic energy density gradient turns out to affect the 
final form of the temperature gradient of the fluid parcel. In 
particular, the relation between the gradient of the parcel and 
that of the surrounding fluid, given by Equation (46). If / = 0, 
then we find that complete adiabaticity holds, meaning V e — > 
V at j. However, for the case that / ^ 0, any heat transferred 
away from the parcel will be in the form of magnetic energy 
(dQ = d%). Thus, from Equation (42), 



Q = c P dT e --dP e 
P 



FSv 
P0CX 



which may be rearranged to read 

V e = V ad 



(51) 



(52) 



Substituting this expression for the parcel's temperature 
gradient back into Equation (47), we derive the condition that 
must be met if a fluid parcel is to be stable against convection, 



V ad -/-V ad V x - (1 -/) gV z > V temp . 



(53) 



With a modified stability criterion in hand, LS95 the proceed 
to derive a set of five equations that allow for a solution to 
the temperature gradient, V temp . The equations are developed 
through a detailed consideration of the various energy fluxes 
through the system in their Sections 5.1 - 5.3. The final five 
equations comprising their magnetic mixing-length descrip- 
tion of convection are 



Ftnt — 



4acG T 4 M r , 

3 KPr 2 
4acG T A M r 



rad 



2acT i 

PVconvCp 



3 

" P^'conv 
CO 

1 +yco 2 



K/V 2 Vtemp " 



cpDT 



temp " 



P8v 
P<*X 

•v.)- 



:(1-/)V, 



(54) 
(55) 
(56) 

(57) 



(V 



temp 



(V e -V ad )+/-Va d V z 



(58) 



where we must now take a moment to dissect the various 
pieces. 

The first two equations above describe the total flux if only 
radiation is carrying energy and in the case that convection 
is also present, respectively. Within those two equations, a, 
c, and G are, respectively, the radiation constant, speed of 
light in a vacuum, and the gravitational constant. M r is the 
mass contained within a spherical volume characterized by 
the radius, r, and K is the radiative opacity. Lastly, V lad is the 
radiative temperature gradient. 

Equation (56) is the energy flux transported by convection. 
The quantities, DT and D% are defined the same as Dp in 
Equation (45). Another variable, the convective velocity v conv 
is also introduced and characterizes the velocity of the fluid 
within a convection cell. 



The convective velocity is then defined in Equation (57) 
and contains a single parameter not previously mentioned, i m . 
This parameter is the convective mixing-length, £ m , which is 
further defined as some multiple of the pressure scale height 
(i.e., t m = Note that the mixing-length introduces 

the canonical convective mixing-length parameter, OCmlt, into 
the discussion. 

Last of the five equations of MLT describes how the con- 
vective gradient is connected to the adiabatic gradient. Here, 
CO = Kp£ m and y is set by the geometry of the convecting bub- 
ble. The shape parameter, y is partially what separates the var- 
ious formulations of MLT. Consistent with the standard DSEP 
treatment of convection, and LS95, we set y = 1 /3. 

Combining Equations (54)-(58) yields a solution for the 
convective velocity, which effectively defines the temperature 
gradient. The entire solution, as with the previous derivations, 
may be found in its full glory in LS95. For our purposes, we 
cite yet another set of new variables required to simplify the 
solution, 

2 = 5 (59) 



Yo 



c P p ( 1 + co 2 /3 



2acT 3 



CO 



C = 



v 



1/2. 



= YoC 

1 /2 

V rad - v ad + /- v ad v z + (i - f)-v x ) 



9 co 
A = - 



(3 + co 2 ) 



and, lastly, 



3' : 



(60) 



(61) 



(62) 



(63) 



(64) 



Resulting from the combination of the five magnetic MLT 
equations and the substitution of the variables just defined 
produces an equation that is quartic in y, 



2A 
~V 



y 4 +y 3 + 



2Af C 



a Q 

y- 



(i 



■/)vV z + l 



Vy 2 



Q 



(1-/)V Z = 0. (65) 



We remark that in the equation above, we deviate from LS95 
by the presence of a 1 /Q in the quadratic term. This difference 
appears to be the result of the authors accidentally dropping a 
term in the original derivation. However, as we shall see, this 
factor will become unimportant. 

Finding a solution for y from Equation (65) can easily be 
obtained numerically. Making an educated guess as to the 
solution for y we may use a series of Newton-Raphson cor- 
rections to converge to a proper solution. A convergence tol- 
erance of 10~ 10 is imposed on the correction term to reduce 
propagation of large numerical errors. Modifying the original 
"guess" supplied by LS95 to include the dropped term men- 
tioned above, a good trial solution is 



l + 2AC^vf^-i)(l-/)V, 



a 



(66) 
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One may notice that in the above trial solution (and most other 
MLT equations), the precise value of the free parameter / has 
the ability to drastically simplify the expression. 

3.4. The Parameter f and the Frozen Flux Condition 

In Section 3.1 we specified that the plasma under consid- 
eration was perfectly conducting and, thus, had zero resis- 
tivity. One consequence of assuming an ideal MHD plasma 
is that magnetic field lines become physical objects that are 
transported by the plasma, the so-called frozen flux condition 
(FFC; Alfven 1942). The magnetic flux is 



*(t) 



S(t) 



R(r,t)-dA. 



(67) 



For an ideal plasma, the evolution of <t> in time is dependent 
upon not only the time rate of change of the magnetic field, 
B(r, t), but also on any distortion occurring to the bound- 
ing surface, S(t), as the plasma moves. The net effect is that 
the time rate of change of the magnetic flux is equal to zero. 
Therefore, we must have 



3B 

a7 



= o, 



(68) 



which may be rewritten using Equation (1), the induction 
equation, with r| = 0. This results in the well-known FFC 
condition, 

V x (u x B) = 0. (69) 

The FFC enforces the restriction that, for a spherically sym- 
metric bubble of plasma undergoing isotropic expansion or 
contraction (Kulsrud 2004), 



B 

373 



constant. 



(70) 



Applying this constraint to the magnetic energy gradient def- 
inition of a convecting fluid element (Equation (46)), we are 
able to physically motivate the definition of the parameter / 
which governs the flux of energy between the fluid element 
and the surrounding material. 

Imagine a region in a star where a small bubble begins to 
grow convectively unstable. Initially, the bubble has the same 
properties as the surrounding fluid. It is only because of the 
change in density that other properties begin changing as well. 
The FFC allows us to write the magnetic energy contained 
within a fluid parcel as a function of the magnetic energy of 
the surrounding material. Since a convecting fluid parcel has 
a slight under- or overdensity compared to its surroundings, 
we perturb the element's density 



P, = P,+^ 



(71) 



where it is understood that \ <C p. s . We also drop the subscript 
s hereafter. Assuming, for simplicity, that the bubble expands 
isotropically, the magnetic field strength within a convectively 
unstable bubble is 



2/3 



(72) 



meaning the magnetic energy per mass may be written as 
Bl B 2 



= 87T(pH) = 8^ 1 + p~, 



(73) 



We now have a direct relation between the magnetic energy 
density of the convecting fluid element and the surrounding 
material, 

Xe = XsU + ^j ■ (74) 

Taking the radial, logarithmic derivative, 

'dln%\ ( d\x\y\ 1 / \ \ \dlnt, c/lnp 



dr 



= ( d\n% 
V dr 



3 VP + ^ 



dr 



dr 



(75) 

The first of the bracketed terms goes to zero, as the density 
perturbation is independent of radial location. Using the defi- 
nition of % s (Equation (20)) to expand the derivative, and after 
rearranging the resulting terms, the equation becomes 



(dln%\ _d in (B 2 \ (d\np\ 
V dr ) e dr \%%) \ dr j 



(76) 

By definition, we know that E,/p <C 1, meaning the perturba- 
tion term in the square brackets is negligible. As such, 



(77) 



allowing us to conclude that the FFC implies that / « 1 . 



3.5. Magnetic Field Strength Distribution 

The strength of the perturbations described in the preceding 
sections are determined by the magnitude and spatial gradient 
of %. Mentioned in Section 3.3, was that we deviate from the 
prescription of % proposed by LS95. Instead of defining 



where 



X — Xmax exp 



M D = log 



1 ( M D -M Dc 

2 V o 



10 



1 



(78) 



(79) 



we opt to directly prescribe the radial magnetic field profile. 
Approaching the problem in this manner, however, introduces 
the difficulty of selecting a particular radial profile, and with- 
out any real confidence of the radial profile inside stars, we 
are left to our own devices. 

One of the simplest profiles to select is that of a dipole con- 
figuration, where the field strength drops off as r~ 3 from the 
magnetic field source location. This is illustrated in Figure 1 . 
The radial profile may then be prescribed as 



B(R) = B(R tach ) 



(R/R, 



tach) 



R > Rtzch 
R < ^tach 



(80) 



with the peak magnetic field strength defined to occur at the 
radius R ta ch- The radial location described by ^ t ach is the lo- 
cation of the stellar tachocline, an interface between the con- 
vective envelope and radiative core. This interface region 
is thought to be characterized as a strong shear layer where 
the differentially rotating convection zone meets the radia- 
tive core rotating as a solid body. Theory suggests that the 
tachocline is the source location for the standard mean-field 
stellar dynamo (i.e., the oc-co dynamo; Parker 1975). 

Since DSEP monitors the shell location of the boundary 
to the convection zone, the tachocline appeared to be a rea- 
sonable location, both theoretically and numerically, to base 
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Figure 1. Magnetic field strength profile for a 1.0 M,., star with a 5.0 kG 
photospheric magnetic field strength (maroon, solid). The green dash-dotted 
line indicates the location of the stellar tachocline, the interface between the 
radiative and convective regions. The plot is meant only to illustrate the field 
strength profile. A small gap in the field strength profile is barely perceptible 
near the surface of the star. This artifact is due to the separation of the stellar 
interior and envelope integration regimes in the code. 

the scaling of the magnetic field strength. However, defining 
the magnetic field strength at the tachocline (B(R tac \ i )) is re- 
quired. In an effort to allow for direct comparisons between 
field strengths input into the code and observed magnetic field 
strengths, the field strength at the tachocline is anchored to the 
photospheric (surface) magnetic field strength, 



B(^tach) 



. "tach / 



(81) 



where B sur f = B(R*) is introduced as a new free parameter. 
The advantage of B sur f as a free parameter is that it has poten- 
tial to be constrained observationally. 

Fully convective stars do not possess a tachocline. How- 
ever, a dynamo mechanism still has the potential to 
drive strong magnetic fields through an a 2 mechanism 
(Chabrier & Kiiker 2006). Full three-dimensional MHD mod- 
eling suggests that, in this case, the magnetic field strength 
peaks at about 0.15 R* (Browning 2008). Unfortunately, 
the adopted micro-physics were solar-like and may not be 
entirely suitable for fully convective M-dwarfs. Regardless 
of these shortcomings, Browning's investigation provides the 
only estimate, to date, for the location of the peak magnetic 
field strength in fully convective stars. As such, we adopt 
^tach = 0.15/?* as the dynamo source location in our models 
of fully convective stars. 

3.6. Numerical Implementation 

Although we have laid out the mathematical construction of 
the magnetic perturbation, we have yet to illuminate precisely 
how the perturbation is treated numerically. When a magnetic 
model is first executed, the user provides a surface magnetic 
field strength, the geometry parameter y, and the age at which 
the magnetic perturbation will occur. The model proceeds to 
evolve the same as a standard model until the initial perturba- 
tion age is reached. 

Once the perturbation age is reached, the magnetic field 
strength profile is prescribed based on the assumed photo- 
spheric field strength and the location of the tachocline, as in 
Figure 1 . The magnetic energy and magnetic pressure are then 



computed for each of the model's mass shells. Here, the total 
pressure associated with each mass shell is also perturbed. 

Following the introduction of the perturbation, the code 
must recompute the structure of the stellar envelope, which 
is separate from the stellar interior integration. The envelope 
comprises the outer 2%-3%, by radius, of the stellar model. 
Surface boundary conditions are determined prior to the enve- 
lope integration by interpolating within the PHOENIX model 
atmosphere tables using logg and T e ff. This interpolation re- 
turns Pijas at the surface of the star, and defines the start of the 
inward integration. The magnetic perturbation is then explic- 
itly included in the calculation of the analytic EOS. 

This leads into the convection routines, where the non- 
standard stability criterion in Equation (53) is evaluated and 
judged. Either the equations of magnetic MET are solved, 
or the radiative gradient is selected. The envelope integra- 
tion scheme proceeds until it reaches a pressure commensu- 
rate with the pressure for the interior regime. 

From the newly calculated envelope, the interior integra- 
tion begins using a Henyey integration scheme (Henyey et al. 
1964) with the magnetic perturbation implemented. The EOS 
and convection routines are evaluated as in the envelope. 
Once a final solution is converged upon, the code iterates in 
time and the process is repeated. For each temporal iteration, 
the magnetic field profile is adjusted to adapt to changes in the 
location of the tachocline and changes in the number of mass 
shells. 

4. INITIAL TESTING 

In Section 3, we outlined the formulation and implemen- 
tation of a magnetic perturbation within the framework of 
DSEP. With the perturbation implemented, it was crucial to 
perform a series of numerical tests and common-sense checks 
to validate that the code was operating properly. 

The four key numerical tests were to ensure that: 

1. The relative change in radius between magnetic mod- 
els of monotonically increasing photospheric magnetic 
field strength must also be monotonically increasing 
with respect to a non-magnetic model. 

2. All model adjustments after the initial perturbation 
must be continuous and smooth. 

3. The final perturbed model properties must be indepen- 
dent of the evolutionary stage at which the perturbation 
is made. 

4. The resulting perturbed model must be consistent, re- 
gardless of the number of time steps taken after the per- 
turbation. 

All of these tests were performed to confirm that the code 
was producing consistent results and that it was doing so in a 
numerically stable manner (i.e., no wild fluctuations). 

Figure 2 demonstrates that all model adjustments to a mag- 
netic perturbation satisfy each of the four criteria we required 
for numerical validation. Panel (a) demonstrates that the ra- 
dius monotonically increases as the surface field strength ap- 
plied monotonically increases. Changes are observed to be 
smooth, as seen in panel (b) and are independent of the num- 
ber (or size) of the evolutionary time steps taken (panel (d)). 
Finally, the plot in panel (c) indicates that the relative change 
to the model asymptotes to the same value, regardless of the 
evolutionary phase at which the perturbation is applied. 
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Figure 2. Tests of numerical stability for a 1 XI . magnetic model with y = 2. (a) Influence of various magnetic field strengths, (b) Evidence of a smooth 
perturbation at a large magnetic field strength with a relatively large radius change, (c) Consistency among models with the perturbation turned on at various 
evolutionary stages, (d) Demonstrating the insensitivity of the perturbation to the number of time steps after the perturbation is enabled. 



Beyond testing for numerical stability, we must be as- 
sured that the code produces results consistent with reality. 
Typically, a comparison with previous studies would be uti- 
lized. However, the only such examples computed for evo- 
lutionary timescales are for CM Dra (Chabrier et al. 2007; 
MacDonald & Mullan 2012). The stellar mass regime oc- 
cupied by CM Dra would require the implementation of 
FreeEOS, a task reserved for a future investigation. With the 
analytical EOS, it would appear there are no models generated 
with which to compare. Even LS95 operate over timescales 
on the order of a solar-cycle and not evolutionary times. 

In the absence of previous results with which to compare, 
we opted to impose a weak magnetic field (5 G) on our solar- 
calibrated model. Seeing as the properties of the Sun do not 
require a magnetic field in order to produce an adequate solar 
model (Bahcall et al. 1997), we would expect the solar prop- 
erties to remain almost entirely unaffected. As expected, the 
magnetically perturbed model still meets our requirements for 
it to be considered solar-calibrated (see Section 2). The model 
radius changes by 3 parts in 10 5 while effective temperature 
changes by 4 parts in 10 6 given the presence of the weak field. 

5. CASE STUDY: EF AQUARII 
5.1. Standard Models 

Standard, non-magnetic mass tracks with solar metallicity 
were computed for both EF Aqr A and B with masses of 
1.24 Mq and 0.95 M@, respectively. Additional mass tracks 
were also generated for a scaled-solar metallicity of +0. 1 dex. 
The two components were unable to be fit with a coeval age, 
regardless of the adopted metallicity. This is consistent with 
the conclusions of Vos et al. (2012). Two fitting methods were 
performed to this end. 

We first attempted to fit both components on an age-radius 



plane, which is equivalent to fitting on the mass-radius (M- 
R) plane. This is illustrated by the solid lines in Figure 3(a). 
The primary star evolves much more rapidly than secondary, 
owing to the rather large mass difference. Thus, the model 
radius of the primary inflates to the observed radius at an age 
of 2.0 Gyr. The radius of the primary then quickly exceeds the 
observational bounds within about 0. 1 Gyr. The observational 
bounds on the primary radius places tight constraints on the 
allowed age our models predict for the system. However, our 
model of the secondary does not reach the observed radius 
until an age of 6.3 Gyr, an age difference of 4.3 Gyr between 
the two components. This 4.3 Gyr difference is consistently 
present in our models, including when mass and composition 
are allowed to vary within the observed limits (not shown). 

Assuming the age of the system was predicted accurately 
using models of EF Aqr A in the M-R diagram, we find that 
the model radius of EF Aqr B underpredicts the observed ra- 
dius by 11%. Again, consistent with the findings of Vos et al. 
(2012), who found a radius discrepancy of 9%. Such a dis- 
agreement is also broadly consistent with results from other 
studies of active EBs (Ribas 2006; Torres et al. 2010). 

A second approach was to fit the system on an HR diagram 
using the individual mass tracks. Figure 3(b) demonstrates 
that the standard model tracks do not match the observed r e g— 
luminosity of either star, despite fitting the stars individually 
in the M-R plane. Our models predict temperatures that are 
250 K (4%) and 430 K (8%) hotter than the observations for 
the primary and secondary, respectively. 

The noticeable temperature disagreements in both stars may 
be the result of two possible effects. On the one hand, we 
might assume that the age predicted from the M-R plane is 
correct and that there only exists a discrepancy with the ef- 
fective temperatures. This implies that either the effective 
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Figure 3. Individual stellar mass tracks representing EF Aqr A and B (labeled). Non-magnetic mass tracks are shown as a red, solid line while y = 2 and y = 4/3 
magnetic tracks are indicated by the blue, long-dash and light-blue, short-dash lines, respectively. The corresponding photospheric magnetic field strengths for 
EF Aqr A are 1.6 kG (y = 2) and 2.6 kG (y = 4/3). For EF Aqr B they are 3.2 kG (y = 2) and 5.5 kG (y = 4/3). (a) Age-radius diagram with the observed radii 
marked as gray regions. Vertical dotted lines define the age constraint imposed by EF Aqr A, suggesting the system has an age of 1.3 ±0.2 Gyr. (b) HR diagram 
with the observed data marked as black points. 



temperature from the models or the observations is incorrect. 
However, since the stars are quite similar to the Sun, it is more 
likely the case that our standard models underpredict the ra- 
dius of the primary as well, driving up the model-derived ef- 
fective temperature. This particular scenario is supported by 
Vos et al. (2012) who found both components display obvi- 
ous Ca II emission that is likely the result of each star having 
a magnetically heated chromosphere. 

The age inferred from the M-R diagram would then be 
older than the true age of the system. Unfortunately, this sce- 
nario further complicates the situation regarding EF Aqr B. 
If the age of the system is younger than inferred from stan- 
dard models of the primary, then the radius discrepancy for 
the secondary becomes larger than originally quoted. 

5.2. Magnetic Models 

Following the results discussed above, magnetic models 
were computed for both EF Aqr components using a scaled- 
solar heavy element composition. Several models with a mass 
of 1.24 Mq were generated with various surface magnetic 
field strengths. We then found the model with the weakest 
field strength required to produce the observed radius and r e ff 
of EF Aqr A. With y = 2, we had to prescribe a photospheric 
magnetic field strength of 1.6 kG, while with y = 4/3, a more 
intense 2.6 kG field was necessary. The magnetic model 
tracks are displayed in Figures 3(a) and (b) as blue dashed 
lines. The magnetic models of the primary suggest a younger 
age of 1 .3 ± 0.2 Gyr for the EF Aqr system, as opposed to the 
2.0 Gyr age determined from standard models. This younger 
age is consistent with the 1 .5 ± 0.2 Gyr age derived for the pri- 
mary by Vos et al. (2012) after fine-tuning the mixing length. 

We next had to select a magnetic field strength that would 
allow a 0.95 Mq model to have a radius and T e s compatible 
with EF Aqr B at 1.3 Gyr, if finding that unique combination 
was possible. Surface magnetic field strengths of 3.2 kG and 
5.5 kG were able to produce the required stellar parameters, 
with y = 2 and 4/3, respectively, at an age of 1.35 Gyr. In both 
cases, the models were able to reproduce the stellar radius and 
r e ff within the quoted la uncertainties. Figures 3(a) and (b) 



demonstrate that the magnetic models do indeed match both 
component radii and r c ffS at a common age. 

Structurally, the addition of a magnetic perturbation within 
the models reduces the radial extent of the surface convection 
zone. For both stars in EF Aqr, we find the magnetic models 
that are sufficient to correct the observed discrepancies have 
surface convection zones that are 4% smaller than those in the 
standard models, at the same age. The reduction of the surface 
convection zone can be attributed to the modified stability cri- 
terion as well as modified convective velocities. While only 
speculation, we attribute the equality of the percent reduction 
of convection zone sizes to coincidence. 

6. DISCUSSION 

6.1. Field Strengths 

The implementation of a magnetic perturbation within stel- 
lar evolution models is quite capable of reconciling predicted 
model fundamental stellar properties with those determined 
observationally, at least for EF Aqr. While it seems plausible 
that magnetic fields may suppress thermal convection inside 
solar-type stars, how are we to be sure that magnetic fields 
may be reasonably invoked for this particular system? Even 
if invoking magnetic fields is rational, are the field strengths 
required by the models realistic? 

Addressing the first question, we showed in Section 1 
that naturally inefficient convection, as described by the 
Bonaca et al. (2012) calibration, was unable to account for 
the small values of OCmlt required to mitigate the observed 
model-observation disagreements. But, the inability of natu- 
rally inefficient convection to provide a solution does not pos- 
itively identify magnetic fields as the root cause. However, 
there is additional evidence that invoking magnetic fields is 
reasonable. 

High-resolution spectroscopy of the Ca II H and K lines 
for both stars in EF Aqr reveals incredibly strong emission 
cores superposed on the absorption troughs (Vos et al. 2012). 
A search of the ROSAT Bright Source Catalogue (Voges et al. 
1999) also shows that EF Aqr is a strong X-ray emitter. Cou- 
pling these observations with high projected rotational veloc- 
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ities extracted from line broadening measurements, suggests 
that there is the potential for a strong dynamo mechanism to 
be operating. This evidence is only circumstantial, but does 
provide tantalizing clues. 

For the sake of argument, let us assume that the stars are sig- 
nificantly magnetically active. It would then be worthwhile to 
compare the strength of the magnetic field for each EF Aqr 
component to those values required by the models. Unfortu- 
nately, no direct magnetic field strength estimates of EF Aqr 
are available, forcing us to base our analysis on indirect mag- 
netic field strength estimates. 

A natural first step would be to compare the EF Aqr com- 
ponents to other known solar-type stars, such as the Sun 
and a Cen A and B (see Table 2). The Sun's mean pho- 
tospheric magnetic field strength is between 0.1 G and 1 
G (Babcock & Babcock 1955; Babcock 1959; Demidov et al. 
2002) with local patches of very intense fields (i.e., sunspots) 
with strengths on the order of 2-3 kG (Hale 1908). Similarly, 
the average longitudinal field strength of a Cen A was deter- 
mined to be less than 0.2 G, after a null detection of a Stokes 
V polarization signature (Kochukhov et al. 201 1). 

The field strengths required by the models of the EF Aqr 
components therefore suggest that the stars are pervaded by 
magnetic fields that typically characterize the intense regions 
of sunspots. This at first appears detrimental to the validity of 
the models. However, studies of the active and quiet Sun, par- 
ticularly sunspot regions, has led to multiple scaling relations 
allowing for an indirect determination of stellar photospheric 
magnetic field strengths. 

One of these scaling relations was observed to exist be- 
tween the X-ray luminosity of an active region and its total un- 
signed magnetic flux (Fisher et al. 1998; Pevtsov et al. 2003). 
The two observables were found to exhibit a power-law rela- 
tion, 

L x <x®P, (82) 

where the power-law index, p, was determined to be 
1.19±0.04 by Fisher et al. (1998). The index was later re- 
vised by Pevtsov et al. (2003) using a more diverse data set, 
including both solar and extrasolar sources. 7 Their revised 
analysis decreased the index to p = 1.15. The magnetic flux 
is defined in the usual manner (Equation (67)), 

4>= [b-cIA = [ \B z \dA. (83) 
Js Js 

with \B Z \ represents the vertical magnetic field strength. 
Therefore, if we are able to determine the X-ray luminosity 
of the EF Aqr components, it is possible to place a lower limit 
on the magnetic field strength at the surface of the two com- 
ponents. 

The system has a confirmed X-ray counter-part in the 
ROSAT All-Sky Survey Bright Source Catalogue (Voges et al. 
1999). The X-ray count rate was converted to an X-ray flux 
according to the formula derived by Schmitt et al. (1995), 

F v = (5.30HR + 8.31)xl0" 12 X cr , (84) 

where HR is the X-ray hardness ratio, 8 X cr is the X-ray count 
rate, and F x is the X-ray flux. Finally, the X-ray flux was 

7 The total unsigned magnetic flux of the stellar sources was obtained using 
direct observational techniques (Saar 1996). 

8 There are typically two hardness ratios listed in the ROSAT catalog, HR1 
and HR2. The Schmitt et al. (1995) formula requires the use of HR1. 



converted to a luminosity using the 172 pc distance quoted by 
Vos et al. (2012). 

The count rate measured by ROSAT was 0.0655 ±0.0 154 
counts s with a hardness ratio of 0.32±0.22. This yields 
an X-ray flux of 6.55 x 10~ 13 erg cm~ 2 . Weighting the con- 
tribution of each star to the total X-ray flux is reliably per- 
formed in one of two ways: by assuming both stars contribute 
equally (valid if for binaries if stars are similar in radius; 
Fleming et al. 1989) or weighting proportional to v ro fSini 
(Pallavicini et al. 1981; Fleming et al. 1989). Given the sim- 
ilarity of the two stars in EF Aqr, the precise weighting does 
not affect the results. If both stars contribute equally to the 
total X-ray flux, then at a distance of 172 pc, the X-ray lu- 
minosity of each component is L x = 1.16 x 10 30 erg s _1 . 
Alternatively, weighting the two stars based on their pro- 
jected rotational velocity gives L X , A = 1.25 x 10 30 ergs -1 and 

L X , B = 1.07 x 10 30 ergs _1 . 

To provide a comparison, the X-ray luminosity of a typical 
solar active region is on the order of 10 27 erg s _1 (Fisher et al. 
1998; Pevtsov et al. 2003). The Sun, on average, has a total 
X-ray luminosity of 10 27 erg s _1 up to nearly 10 28 erg s _1 , 
depending on where in its activity cycle it is located (Ayers 
2009). Similarly, Ayers (2009) monitored the X-ray luminos- 
ity of a Cen and found the primary had an X-ray luminosity 
around half that of the Sun (approximately 10 27 erg s _1 ) and 
the secondary had about twice the X-ray luminosity of the 
Sun (about 10 28 erg s _1 ). Further estimates for the X-ray lu- 
minosity of a Cen A and B comes from ROSAT, which yields 
luminosities between 10 27 erg s and 10 28 erg s _1 for each 
component, consistent with Ayers' analysis. Table 2 provides 
a comparison of how these quantities translate to magnetic 
field strengths. 

Comparisons with the Sun and a Cen show that each com- 
ponent in EF Aqr has an X-ray luminosity 2-3 orders of mag- 
nitude greater than "typical" G and early-K stars. Again, 
while not indicative of causation, the correlation between high 
levels of X-ray emission and magnetic activity strongly sug- 
gests that EF Aqr is incredibly active. Given this information, 
our initial assumption that the stars are active seems valid. 
Therefore, the implementation of a magnetic perturbation in 
stellar models of this system appears warranted. 

The amount of vertical magnetic flux near the surface of 
each star may then be found using the Pevtsov et al. (2003) 
scaling relation. This suggests <t> = 1.39 x 10 26 Mx asso- 
ciated with each component (given equal flux contribution). 
Converting to a magnetic field strength involves dividing the 
unsigned magnetic flux by the total area through which the 
field is penetrating. For our purposes, the area is the entire sur- 
face area of the star. Therefore, we find the vertical magnetic 
field strength for the primary and secondary of EF Aqr to be 
1.3 kG and 2.5 kG, regardless of the adopted flux weighting, 
respectively. Note, this is the vertical magnetic field strength 
and sets a lower limit to the total magnetic field strength. It 
should also be mentioned that the X-ray luminosities calcu- 
lated for EF Aqr A and B are near the edge of the data sample 
utilized by Pevtsov et al., although no extrapolation of the re- 
lation was required. 

Further estimates of the magnetic field strengths may be 
found by applying a scaling relation using the Ca II K line 
core emission (Schrijver et al. 1989). The scaling relation was 
developed by correlating Ca II K line core emission and the 
magnetic flux density of solar active regions and their sur- 



Magnetic Stellar Models of EF Aquarii 



13 



Table 2 

Estimated Magnetic Field Strengths (in G) 



Star 


Direct 


X-Ray 


Call 


DSEP 


Sun 


0.1 - 1.0 


5-20 






a Cen A 


<0.2 


~3 






a Cen B 




~49 






EF Aqr A 




1300 


830 


1600 - 2600 


EF Aqr B 




2500 


3300 


3200 - 5500 



roundings. The relation between the two was found to be 



— — — — = 0.008 (B/) 0,6 (85) 

where I c is the intensity of the emission line core, I w is the 
intensity of the line wing. While the relation was derived for 
small, local active regions, Schrijver et al. (1989) suggest that 
there is no reason to believe that the relation would not hold 
for hemispherical averages of solar-type stars. 

Using the spectra provided by Vos et al. (2012) for the Ca II 
K core emission lines from each EF Aqr component, we were 
able to estimate the magnetic field strength of each compo- 
nent. Spectra for the primary indicate that the average mag- 
netic field strength, (Bf), is equal to 830 G. Similarly, for 
EF Aqr B, (Bf) = 3.3 kG. The values quoted above are derived 
from rough approximations of the core and wing intensities. 
However, we do not foresee the values of (Bf) changing rad- 
ically with more precise line intensity measurements. We do 
caution that the results for EF Aqr B require an extrapolation 
of the Ca II relation and the data for EF Aqr A place it near 
the edge of the derived relation where only a few data points 
exist. 

Based on the scaling relations for X-ray emission and 
Ca II K line core emission, the magnetic field strength for 
the primary and secondary is seen to be approximately 1 kG 
and 3 kG, respectively. The magnetic field strengths required 
by the models are therefore within a factor of two of the pre- 
dicted field strengths, regardless of the adopted y value. Since 
we expect the X-ray emission prediction to be a lower limit to 
the full magnetic field strength, this is extremely encouraging. 
The models do not require abnormally large field strengths to 
reconcile the model properties with those from observations, 
particularly when a y value of 2 is adopted. 

6.2. Implications 

The introduction of self-consistent magnetic stellar evo- 
lution models has multiple applications, ranging from stud- 
ies of exoplanet host stars (Torres 2007; Charbonneau et al. 
2009; Muirheadet al. 2012) to investigations of cataclysmic 
variable (CV) donor stars (Knigge et al. 2011), as well 
as to studies attempting to probe the stellar initial mass 
function of young clusters, where stars are typically very 
magnetically active (Johns-Krull 2007; Jackson et al. 2009; 
Yang & Johns-Krull 2011). Although, most obvious, are the 
implications for studies of low-mass eclipsing binary sys- 
tems (see, e.g., Torres et al. 2010; Parsons et al. 2012, and 
references therein). Low-mass stellar evolution models have 
been highly criticized for being unable to predict the radii 
and effective temperatures of DEB stars. Models incorporat- 
ing magnetic effects open the door to probing the underlying 
cause of the model-observation disagreements and providing 
semi-empirical corrections to models. 

Magnetic fields have long been theorized as the culprit, but 



previous generations of models have only treated magnetic 
fields in an ad hoc manner. Comparing the results of these 
methods with the one presented in this work, both in terms of 
surface parameters and the underlying interior structure, will 
provide an interesting test of their validity. Ultimately, the ad 
hoc models disagree on the dominant physical mechanism un- 
derlying the observed discrepancies. The availability of self- 
consistent magnetic models should help to settle the debate as 
to which mechanism (suppressed convection or starspots) is 
most at work. 

For EF Aqr, the models suggest that magnetic suppression 
of thermal convection is sufficient to reconcile the models 
with the observations. Since stars with small convective en- 
velopes, such as those discussed in this work, are more sen- 
sitive to adjustments of the convective properties, it is not 
wholly surprising that suppressing convection is sufficient to 
explain the observations. Whether this mechanism will be ad- 
equate for stars near the fully convective boundary has yet to 
be seen. Future work modeling the lowest-mass DEB sys- 
tems will clarify this ambiguity. Regardless, this may suggest 
why the largest radius deviations are predominantly observed 
at higher masses (Feiden & Chaboyer 2012), with the notable 
exception of CM Dra (Terrien et al. 2012). 

The nature of our models allows for independent verifica- 
tion of the magnetic field strengths required as input. While 
the indirect estimates provided by X-ray emission and Ca II K 
emission are encouraging, confirmation of these results using 
high resolution Zeeman spectroscopy, spectropolarimetry, or 
Zeeman-Doppler imaging (ZDI) is preferred. Unfortunately, 
these observations are difficult for fast rotators, such as those 
that comprise most DEB systems. They are also difficult for 
distant systems, where the short integration time required by 
ZDI inhibits the ability of acquiring measurements with suffi- 
cient signal-to-noise (see the reviews by Donati & Landstreet 
2009; Reiners 2012). Once a magnetic field is detected, there 
exists the question of whether the observed strength is indica- 
tive of the total magnetic field strength. Field strengths de- 
rived for stars with spectral-type K and M using Stokes V 
observations appear to yield only around 10% of the total 
magnetic field strength compared to observations in Stokes 
/ (Reiners & Basri 2009). This is a consequence of the fact 
that regions of opposite polarity tend to cancel out in Stokes 
V, making it most sensitive to the large-scale component, not 
the small-scale fields thought to pervade low-mass stars. How 
to accurately account for this when testing the models is not 
fully clear and will require investigation. As more stars across 
all spectral types are observed in both Stokes V and /, a more 
coherent picture is sure to develop. 

Magnetic models may also be useful for transiting ex- 
oplanet surveys, particularly those focused on M-dwarfs 
(e.g., MEarth transit survey; Nutzman & Charbonneau 2008; 
Irwin et al. 2009). One of the largest uncertainties in deriv- 
ing the properties of a transiting planet is the radius of the 
host star. The lack of reliability involved in predicting low- 
mass stellar radii from evolution models has deterred the use 
of models as predictors of exoplanet host-star radii (Torres 
2007; Charbonneau et al. 2009). Shoring up these deficien- 
cies may lead to more accurate predictions of host star radii 
from stellar models, circumventing, for a time, the need for 
lengthy and costly observations. This would be most useful in 
identifying interesting follow-up targets by providing a better 
estimate of the habitable zone (Muirhead et al. 2012). 

There are certainly caveats with models, as other large un- 
certainties exist in predicting the properties of a single star 
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from stellar evolution mass tracks (Basu et al. 2012). How- 
ever, work is being performed to alleviate some of these 
uncertainties by calibrating models to asteroseismic data 
(Bonacaetal. 2012). Low-mass stars are also less sensitive 
to the input parameters of stellar models than their solar-type 
counterparts, reducing the associated uncertainties. Stellar 
evolution models may therefore provide a fast and reliable es- 
timate of the host star properties, depending on the required 
level of precision. 

Since most M-dwarfs being surveyed are nearby, there is 
a good chance that they may have an X-ray counterpart in 
either the Bright Source or Faint Source Catalogue from the 
ROSAT All-Sky Survey (Voges et al. 1999, 2000). As was 
demonstrated in Section 6. 1, magnetic field strengths required 
by the models are within about a factor of two (or better, if y = 
2 is adopted) of those predicted by the X-ray scaling relation 
of Pevtsov et al. (2003). This will allow an intelligent choice 
of the magnetic field strength used as an input for the models, 
thus producing more reliable results from stellar models. 

All told, the introduction of a self-consistent set of magnetic 
stellar evolution models provides the potential for models to 
be used with greater reliability in a wide range of applica- 
tions. There still exist several challenges that require attention 
(Boyajian et al. 2012), but this is a first step in addressing key 
issues that have been raised in the past two decades. 
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